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LONG-TERM  GOALS 

The  long-term  technical  goal  is  to  design,  develop  and  test  the  next  generation,  primitive  equation 
ocean  model  for  high-resolution  scientific  (ROMS:  Regional  Ocean  Modeling  System)  and  operational 
(TOMS:  Terrain-following  Ocean  Modeling  System)  applications.  This  project  will  improve  the  ocean 
modeling  capabilities  of  the  U.S.  Navy  for  relocatable,  coastal,  coupled  atmosphere-ocean  forecasting 
applications.  It  will  also  benefit  the  ocean  modeling  community  at  large  by  providing  the  current  state- 
of-the-art  knowledge  in  physics,  numerical  schemes,  and  computational  technology. 

OBJECTIVES 

The  main  objective  is  to  produce  a  tested  expert  ocean-modeling  framework  for  scientific  and 
operational  applications  over  a  wide  range  of  spatial  (coastal  to  basin)  and  temporal  (days  to  seasons) 
scales.  The  primary  focus  is  to  implement  the  most  robust  set  of  options  and  algorithms  for  relocatable 
coastal  forecasting  systems  nested  within  basin-scale  operational  models  for  the  Navy.  The  system 
includes  some  of  the  analysis  and  prediction  tools  that  are  available  in  Numerical  Weather  Prediction 
(NWP),  such  as:  4-dimensional  variational  data  assimilation  (4D-Var),  ensemble  prediction,  adaptive 
sampling,  and  circulation  stability  and  sensitivity  analysis. 

APPROACH 

The  structure  of  TOMS  is  based  on  ROMS  because  of  its  accurate  and  efficient  numerical  algorithms, 
tangent  linear  and  adjoint  models,  variational  data  assimilation,  modular  coding  and  explicit  parallel 
structure  conformal  to  modem  computer  architectures  (both  cache-coherent  shared-memory  and 
distributed  cluster  technologies).  Currently,  both  ROMS  and  TOMS  are  identical  and  continue 
improving  and  evolving.  ROMS  remains  as  the  scientific  community  model  while  TOMS  becomes  the 
operational  community  model. 

ROMS/TOMS  is  a  three-dimensional,  free-surface,  terrain-following  ocean  model  that  solves  the 
Reynolds-averaged  Navier-Stokes  equations  using  the  hydrostatic  vertical  momentum  balance  and 
Boussinesq  approximation  (Haidvogel  et  al.  2000,  2008;  Shchepetkin  and  McWilliams,  2005,  2009). 
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The  governing  dynamical  equations  are  discretized  on  a  vertical  coordinate  that  depends  on  the  local 
water  depth.  The  horizontal  coordinates  are  orthogonal  and  curvilinear  allowing  Cartesian,  spherical, 
and  polar  spatial  discretization  on  an  Arakawa  C-grid.  Its  dynamical  kernel  includes  accurate  and 
efficient  algorithms  for  time-stepping,  advection,  pressure  gradient  (Shchepetkin  and  McWilliams 
2003,  2005),  several  subgridscale  parameterizations  (Durski  et  al.,  2004;  Warner  et  al.,  2005)  to 
represent  small-scale  turbulent  processes  at  the  dissipation  level,  and  various  bottom  boundary  layer 
formulations  to  determine  the  stress  exerted  on  the  flow  by  the  bottom. 

Several  adjoint-based  algorithms  exist  to  explore  the  factors  that  limit  the  predictability  of  the 
circulation  in  regional  applications  for  a  variety  of  dynamical  regimes  (Moore  et  al.,  2004,  2009). 

These  algorithms  use  the  ideas  of  Generalized  Stability  Theory  (GST)  in  order  to  identify  the  most 
unstable  directions  of  state-space  in  which  errors  and  uncertainties  are  likely  to  grow.  The  resulting 
singular  vectors  can  be  used  to  construct  ensembles  of  forecasts  by  perturbing  initial  and  boundary 
conditions  (optimal  perturbations)  and/or  surface  forcing  (stochastic  optimals).  Perturbing  the  system 
along  the  most  unstable  directions  to  the  state-space  yields  infonnation  about  the  first  (ensemble  mean) 
and  second  (ensemble  spread)  moments  of  the  probability  density  function.  Given  an  appropriate 
forecast  skill  measure,  the  circulation  is  predictable  if  low  spread  and  unpredictable  if  large  spread. 

ROMS/TOMS  uniquely  supports  three  different  4D-Var  data  assimilation  methodologies  (Moore  et  al., 
2011a,  b):  a  primal  fonn  of  the  incremental  strong  constraint  4D-Var  (I4D-Var),  a  strong/weak 
constraint  dual  form  of  4D-Var  based  on  the  Physical-space  Statistical  Analysis  System  (4D-PSAS), 
and  a  strong/weak  constraint  dual  form  of  4D-Var  based  on  the  indirect  representer  method  (R4D-Var). 
In  the  dual  fonnulations,  the  search  for  the  best  ocean  circulation  estimate  is  in  the  subspace  spanned 
only  by  the  observations,  as  opposed  to  the  full  space  spanned  by  the  model  as  in  the  primal 
fonnulation.  Although  the  primal  and  dual  formulations  yield  identical  estimates  of  the  ocean 
circulation  for  the  same  a  priori  assumptions,  there  are  practical  advantages  and  disadvantages  to  both 
approaches  (Moore  et  al.,  2011a,  b,  c).  To  our  knowledge,  ROMS/TOMS  is  the  only  open-source, 
ocean  community-modeling  framework  supporting  all  these  variational  data  assimilation  methods  and 
other  sophisticated  adjoint-based  algorithms. 

There  are  several  biogeochemical  models  available  in  ROMS.  In  order  of  increasing  ecological 
complexity  these  include  three  NPZD-type  models  (Franks  et  al.,  1986;  Powell  et  al.,  2006;  Fiechter  et 
al.,  2009),  a  nitrogen-based  ecosystem  model  (Fennel  et  al.,  2006,  2008),  a  Nemuro-type  lower  level 
ecosystem  model  (Kishi  et  al.,  2007),  and  a  bio-optical  model  (Bissett  et  al.,  1999). 

ROMS  includes  a  sediment-transport  model  with  an  unlimited  number  of  user-defined  cohesive  (mud) 
and  non-cohesive  (sand)  sediment  classes  (Warner  et  al.,  2008).  Each  class  has  attributes  of  grain 
diameter,  density,  settling  velocity,  critical  stress  threshold  for  erosion,  and  erodibility  constant.  A 
multi-level  bed  framework  tracks  the  distribution  of  every  size  class  in  each  layer  and  stores  bulk 
properties  including  layer  thickness,  porosity,  and  mass,  allowing  the  computation  of  bed  morphology 
and  stratigraphy.  Also  tracked  are  bed-surface  properties  like  active-layer  thickness,  ripple  geometry, 
and  bed  roughness.  Bedload  transport  is  calculated  for  mobile  sediment  classes  in  the  top  layer. 

ROMS  is  a  very  modem  and  modular  code  written  on  F90/F95.  It  uses  C-preprocessing  to  activate  the 
various  physical  and  numerical  options.  The  parallel  framework  is  coarse-grained  with  both  shared- 
memory  (OpenMP)  and  distributed-memory  (MPI)  paradigms  coexisting  in  the  same  code.  Because  of 
its  construction,  the  parallelization  of  the  adjoint  is  only  available  for  MPI.  Several  coding  standards 
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have  been  established  to  facilitate  model  readability,  maintenance,  and  portability.  All  the  state  model 
variables  are  dynamically  allocated  and  passed  as  arguments  to  the  computational  routines  via 
dereferenced  pointer  structures.  All  private  arrays  are  automatic;  their  size  is  determined  when  the 
procedure  is  entered.  This  code  structure  facilitates  computations  over  nested  grids. 

WORK  COMPLETED 

A  complete  overhaul  of  the  ROMS/TOMS  4D-Var  data  assimilation  capabilities  was  done  in  FY09.  In 
FY10,  we  concentrated  on  developing  and  testing  additional  algorithms  for  4D-Var  observations 
impact,  observations  sensitivity,  array  modes,  and  posterior  error  covariance  analysis.  Some  of  these 
developments  required  the  adjoint  of  the  4D-Var  data  assimilation  system,  denoted  as  (4DVar)T. 
Currently,  (4DVar)T  is  only  available  for  the  dual  formulation  algorithms  (4D-PSAS  and  R4D-Var) 
because  the  adjoint  of  the  Lanczos-based,  conjugate  gradient  algorithm  used  in  the  minimization  is 
much  simpler.  This  capability  is  challenging  in  I4D-Var  due  to  the  complexity  of  the  I/O  intensive, 
Lanczos-based,  conjugate  gradient  and  preconditioning  algorithms  that  are  used  when  the  minimization 
is  carried  out  in  model  space.  We  are  planning  to  build  (I4D-Var)T  in  the  future. 

The  observation  impact  algorithm  can  be  used  to  quantify  the  contribution  of  each  observation  during  a 
4D-Var  analysis  to  a  specified  aspect  (scalar  function,  say  7)  of  the  ocean  circulation  (Langland  and 
Baker,  2004;  Gelaro  and  Zhu,  2009;  Tremolet,  2008;  Moore  et  al.,  2011c).  That  is,  it  identifies  the  part 
of  the  model  space  that  controls  /  and  that  is  activated  by  the  observations.  It  yields  the  actual 
contribution  of  each  observation  to  the  circulation  increment. 

The  observation  sensitivity  is  based  on  (4D-Var)  and  quantifies  the  change  that  would  result  in  the 
circulation  estimate  as  result  of  changes  in  the  observations  or  observation  array  (Tremolet,  2008; 
Moore  et  al,  2011c).  It  is  a  very  useful  tool  for  efficient  generation  of  observation  system  experiments 
(OSEs),  observation  array  design,  and  adaptive  sampling.  It  can  also  be  used  to  predict  the  changes  that 
will  occur  in  I  in  the  event  of  a  platform  failure  or  a  change  in  the  observation  array  (Moore  et  al. , 
2011c). 

The  array  modes  of  the  stabilized  representer  matrix  can  be  used  to  determine  the  most  stable 
component  of  the  circulation  with  respect  to  changes  in  the  innovation  vector  (Bennett,  1985;  Moore  et 
al.,  2011c).  The  array  modes  are  independent  of  the  observation  values  and  depend  only  on  the 
observation  locations,  the  prior  covariances,  and  prior  circulation. 

A  ROMS  training  and  a  numerical  ocean  modeling  workshop  was  held  at  the  Hong  Kong  University  of 
Science  and  Technology  (HKUST),  Honk  Kong,  China,  in  collaboration  with  Dale  Haidvogel 
(Rutgers),  Kate  Hedstrom  (University  of  Alaska),  and  Jianping  Gan  (HKUST)  from  January  5-14, 
2011. 
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RESULTS 


Three  types  of  nesting  capabilities  have  been  designed  in  ROMS:  (i)  refinement  grids  which  provide 
increased  resolution  (3:1,  5:1,  or  7:1)  in  a  specific  region;  (ii)  mosaics  which  connect  several  grids 
along  their  edges,  and  (iii)  composite  grids  which  allow  overlap  regions  of  aligned  and  non-aligned 
grids.  The  mosaic  and  composite  grid  code  infrastructures  are  identical.  The  differences  are 
geometrical  and  primary  based  on  the  alignment  between  adjacent  grids.  All  the  mosaic  grids  are 
exactly  aligned  with  the  adjacent  grid.  In  general,  the  mosaic  grids  are  special  case  of  the  composite 
grids. 

An  example  of  refinement  nested-grids  is  shown  in  Fig.  1  for  the  Monterrey  Canyon  system  with 
increasing  horizontal  resolution.  The  nested  grids  are  imbedded  in  a  regional,  eddy  resolving  grid  of 
the  US  west  coast  (Fig.  la)  with  an  average  resolution  of  7.5km.  The  intermediate  grid  (Fig.  lb)  has  an 
average  resolution  of  2.5km  and  it  is  primarily  used  to  simulate  the  submesoscale  circulation  from 
north  of  Point  Conception  to  Point  Arenas.  The  refinement  ratio  is  1:3.  A  fine-resolution  grid  of 
Monterrey  Bay  is  shown  in  Fig.  lc  with  an  average  resolution  of  500m  and  is  used  to  study  the 
complex  canyon  circulation.  The  grid  lc  has  a  5:1  refinement  ratio  from  lb.  This  type  of  grid  nesting 
has  been  available  to  atmospheric  models  for  several  years  and  more  recently  to  ocean  models. 
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Figure  1:  A  ROMS  nested  grid  approach,  with  increasing  resolution  for  the  California  Current 
System.  The  average  horizontal  resolutions  are:  (a)  7.5km,  (b)  2.5km,  a  nested  ratio  1:3,  and  (c) 
500m,  a  nested  ratio  1:5.  The  grids  mesh  has  been  removed  for  clarity. 


Figure  2  shows  a  nested  application  for  the  US  east  coast  (grid  a)  with  two  imbedded  composite  grids 
for  the  Delaware  Bay  (grid  b)  and  Chesapeake  Bay  (grid  c)  estuaries.  In  addition,  it  has  an  imbedded 
fine  resolution  refinement  grid  d  over  the  Nantucket  Shoals.  This  example  shows  the  unique 
capabilities  of  multiple  grids  nesting  in  ROMS  design.  All  the  grids  can  coexist  in  the  same 
application  and  are  nested  synchronously  at  every  time  step.  This  level  of  complexity  in  nesting  design 
infrastructure  required  a  careful  overhaul  of  ROMS  numerical  kernels  and  extensive  testing  over  the 
past  couple  of  years.  However,  we  have  been  maturing,  testing  and  getting  ready  for  this  development 
since  ROMS  3.0  was  released  in  April  of  2007.  John  Warner  and  I  have  been  working  on  and  off  with 
this  for  several  years.  Indeed,  to  our  knowledge  this  is  one  of  the  most  complex  and  generic  nesting 
approaches  in  an  ocean  model.  Therefore,  the  implementation  and  porting  to  the  community  version  of 
the  code  was  divided  in  three  sequential  development  phases: 
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Figure  2:  A  US  east  coast  example  of  combined  composite  and  refinement  nested-grids:  (a)  parent 
large  scale,  mesoscale  resolution  grid,  (b)  Delaware  Bay  estuary  fine  resolution  composite  grid,  (c) 
Chesapeake  Bay  estuary  system  fine  resolution  composite  grid,  and  (d)  Nantucket  Shoals  fine 
resolution  refinement  grid.  The  grid  meshes  in  (a),  (b),  and  (c)  have  been  removed  and  (d)  has  been 

substantially  sampled  for  clarity. 


•  Phase  I:  the  numerical  kernel  of  ROMS  nonlinear  (NLM),  tangent  linear  (TLM),  representer 
(RPM),  and  adjoint  (ADM)  models  was  modified  to  allow  different  I-  and  J-ranges  in  the  DO- 
loops  to  permit  operations  on  various  nested  grid  classes  ( refinement ,  mosaics,  and  composite) 
and  nesting  layers  ( refinement  and  composite  grid  combinations).  This  facilitates  the 
computation  of  any  horizontal  operator  (advection,  diffusion,  gradient,  etc.)  in  the  nesting 
overlap  regions  and  avoids  the  need  for  cumbersome  lateral  boundary  conditions  on  the  model 
variables  and  their  associated  flux/gradient  values.  The  advantage  of  this  approach  is  that  it  is 
generic  to  any  discrete  horizontal  operator.  Otherwise,  we  have  to  specify  logic  for  each 
advection  scheme  available  in  ROMS.  The  overlap  region  is  an  extended  section  of  the  grid 
that  overlays  an  adjacent  grid.  The  strategy  is  to  compute  the  full  horizontal  operator  at  the 
contact  points  between  nested  grids  instead  of  specifying  boundary  conditions,  as  shown  in  Fig. 
1  of  Warner  et  al.  (2010).  The  stencil  footprint  at  the  contact  edge(s)  of  the  child  grid  is 
expanded  up  to  five  additional  computational  grid-points  to  allow  the  full  evaluation  of  the 
horizontal  operator  at  its  physical  boundary.  The  values  in  the  overlap  region  are  determined 
from  the  parent  grid.  In  refinement  applications,  the  overlap  values  are  computed  by  spatial 
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averaging.  This  is  because  the  parent  grid  is  always  larger  and  coarser  than  the  child  grid  by  a 
factor  of  3,  5,  or  7.  In  mosaic  and  composite  grids,  the  overlap  values  are  computed  by 
interpolation  from  the  parent  grid.  If  child  and  parent  have  coincident  grid  points  in  the 
overlap  region,  the  interpolation  is  not  necessary  and  identical  solutions  are  obtained  when 
compared  to  that  of  one  single  large  and  continuous  grid  (Warner,  et  ah,  2010).  A  coincident 
grid  point  is  one  with  the  same  cell  size,  area,  orientation,  bathymetry,  and  land/sea  masking 
properties.  In  such  cases,  the  coincident  points  are  only  needed  in  the  overlap  region. 

•  Phase  II:  the  lateral  boundary  conditions  C-preprocessing  options  were  eliminated  and 
replaced  with  logical  switches  that  depend  on  the  nested  grid,  if  any.  This  facilitates,  in  a 
generic  way,  the  processing  or  not  of  lateral  boundary  conditions  in  applications  with  nested 
grids.  As  mentioned  above,  the  values  at  the  lateral  boundary  points  are  computed  directly  in 
the  overlap  region  by  the  numerical  kernel.  In  addition,  the  logical  switches  allow  different 
lateral  boundary  conditions  types  between  active  (temperature  and  salinity)  and  passive 
(biology,  sediment,  inert,  etc.)  tracers.  The  lateral  boundary  condition  switches  for  each  state 
variable  and  boundary  edge  are  now  specified  in  ROMS  input  script  file,  ocean. in. 

•  Phase  III:  the  nesting  calls  will  be  added  to  ROMS  main  time-stepping  routines,  main2d  and 
main3d.  The  routine  main2d  is  only  used  in  shallow-water  2D  barotropic  applications 
whereas  main3d  is  used  in  full  3D  baroclinic  applications.  These  main  time-stepping  routines, 
in  turn,  call  several  routines  of  the  ROMS  kernel  to  sub-time  step  (predictor  and  corrector)  each 
tenn  in  the  governing  equations  sequentially:  right-hand-side  terms,  turbulent  mixing,  2D 
momentum,  3D  momentum,  continuity  equation,  and  tracers.  In  Phase  III,  several  routines 
will  be  added  to  process  the  infonnation  that  it  is  required  in  the  overlap  region,  what 
infonnation  needs  to  be  exchanged  from/to  another  grid,  and  when  to  exchange  it.  In  grid 
refinement  applications,  the  information  is  exchanged  at  the  end  of  the  full  time-step  (bottom  of 
main2d  or  main3d).  Contrarily,  in  mosaic  and  composite  grid  applications,  the  infonnation  is 
exchanged  between  each  sub-time  step  call  in  main2d  or  main3d.  That  is,  the  parent  and  the 
mosaic! composite  grids  need  to  sub-time  step  the  2D  momentum  equations  before  any  of  them 
start  solving  and  coupling  the  3D  momentum  equations.  Therefore,  the  governing  equations  are 
solved  and  nested  in  a  synchronous  fashion.  The  concept  of  nesting  layers  is  introduced  to 
allow  applications  with  both  composite  grids  and  refinement  grids,  as  shown  in  Fig.  2.  Here, 
two  nesting  layers  are  required.  In  the  first  nesting  layer,  the  parent  grid  (a)  and  composite 
grids  (b  and  c)  are  synchronously  sub-time  stepped  for  each  governing  equation  tenn.  Then,  in 
the  second  and  last  nesting  layer,  the  refinement  grid  (d)  is  time  stepped  and  the  infonnation 
between  parent  (a)  and  refinement  (d)  grids  are  exchanged  at  the  end  of  the  time  step.  The 
exchange  between  grids  can  be  one-way  or  two-way. 

The  changes  made  to  the  code  during  Phases  I  and  II  were  significant  and  required  several  months  of 
coding  and  over  a  year  to  extensively  test  all  the  nonlinear  and  adjoint-based  algorithms.  All  the 
idealized  and  realistic  test  cases  distributed  with  ROMS  were  updated  and  tested  to  guarantee  identical 
solutions  with  previous  versions  of  ROMS  in  both  serial  and  parallel  (MPI  and  OpenMP) 
computations.  Phase  I  was  released  to  the  community  as  ROMS  3.5  on  April  25,  2011.  Phase  II  was 
released  as  ROMS  3.6  on  September  23,  201 1.  The  final  Phase  III  will  be  released  as  ROMS  4.0  in 
the  near  future,  hopefully  sometime  in  January  2012. 
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A  Graphical  User  Interface  (GUI)  is  currently  under  development  to  help  prepare  the  ROMS  standard 
input  fde,  ocean. in.  The  added  complexity  in  ROMS  nested  grids  makes  this  input  fde  very  confusing 
and  difficult  to  edit  by  hand.  A  sample  of  the  ROMS  GUI  is  shown  in  Fig.  3.  The  GUI  is  written  in 
Python  to  facilitate  cross  platfonn  capabilities,  ease  of  maintenance,  and  code  readability.  This  GUI 
will  be  released  to  the  community  in  the  near  future. 


Figure  3:  ROMS  standard  input  file  python-based  Graphical  User  Interface  (GUI). 
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IMPACT/APPLICATIONS 


This  project  will  provide  the  ocean  modeling  community  with  a  freely  accessible,  well  documented, 
open-source,  terrain-following,  ocean  model  for  regional  nowcasting  and  forecasting  that  includes 
advanced  data  assimilation,  ensemble  prediction,  and  analysis  tools  for  adaptive  sampling  and 
circulation  dynamics,  stability,  and  sensitivity. 

TRANSITIONS 

The  full  transition  of  ROMS/TOMS  to  the  operational  community  is  likely  to  occur  in  the  future. 
However,  the  ROMS/TOMS  algorithms  are  now  available  to  the  developers  and  scientific  and 
operational  communities  through  the  website  http://www.myroms.org/. 

RELATED  PROJECTS 

The  work  reported  here  is  related  to  other  already  funded  ONR  projects  using  ROMS.  In  particular,  the 
PI  (H.  Arango)  closely  collaborates  with  A.  Moore  (adjoint-based  algorithms)  at  University  of 
California,  Santa  Cruz,  A.  Miller  and  B.  Comuelle  (ROMS  adjoint  and  variational  data  assimilation)  at 
Scripps  Institute  of  Oceanography,  E.  Di  Lorenzo  (Southern  California  predictability)  at  Georgia 
Institute  of  Oceanography,  and  J.  Wilkin  (Mid-Atlantic  Bight  variational  data  assimilation)  at  Rutgers 
University. 

The  PI  is  also  supported  by  the  following  grants: 

“Characterization  and  Modeling  of  the  Philippine  Archipelago  Dynamics  Using  the  ROMS  4DVAR 
Data  Assimilation  System”,  grant  number  N00014-04-1-0417,  http://www.myroms.org/philex. 
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